### produces Figure 1
library(gplots)
library(foreign)
library(MASS)
#setwd("e:/info/") 
ciriop<-read.csv("ciriopf1rn.csv",header=F,sep=",")
names(ciriop)
xx<-seq(0:8)

#pdf("e:/tex/a/cwhrf1rn.pdf")
postscript("e:/tex/a/cwhrf1rn.eps")
par(mfrow=c(3,3))


plotCI(ciriop$V1,xx,li=ciriop$V3,ui=ciriop$V4,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 0",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V5,xx,li=ciriop$V7,ui=ciriop$V8,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 1",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V9,xx,li=ciriop$V11,ui=ciriop$V12,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 2",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V13,xx,li=ciriop$V15,ui=ciriop$V16,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 3",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V17,xx,li=ciriop$V19,ui=ciriop$V20,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 4",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V21,xx,li=ciriop$V23,ui=ciriop$V24,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 5",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V25,xx,li=ciriop$V27,ui=ciriop$V28,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 6",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V29,xx,li=ciriop$V31,ui=ciriop$V32,xaxt="n",gap=0,ylab="physical integrity",
xlim=c(-1,1),axes=F,err="x",main="previous value 7",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)

plotCI(ciriop$V33,xx,li=ciriop$V35,ui=ciriop$V36,xaxt="n",gap=0,ylab="physical integrity",
xlim=c(-1,1),axes=F,err="x",main="previous value 8",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)

dev.off()

### produces Figure 2 


ciriop<-read.csv("ciriopf1qn.csv",header=F,sep=",")
names(ciriop)
xx<-seq(0:8)

#pdf("e:/tex/a/cwhrf1qn.pdf")
postscript("e:/tex/a/cwhrf1qn.eps")
par(mfrow=c(3,3))


plotCI(ciriop$V1,xx,li=ciriop$V3,ui=ciriop$V4,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 0",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V5,xx,li=ciriop$V7,ui=ciriop$V8,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 1",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V9,xx,li=ciriop$V11,ui=ciriop$V12,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 2",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V13,xx,li=ciriop$V15,ui=ciriop$V16,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 3",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V17,xx,li=ciriop$V19,ui=ciriop$V20,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 4",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V21,xx,li=ciriop$V23,ui=ciriop$V24,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 5",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V25,xx,li=ciriop$V27,ui=ciriop$V28,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 6",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V29,xx,li=ciriop$V31,ui=ciriop$V32,xaxt="n",gap=0,ylab="physical integrity",
xlim=c(-1,1),axes=F,err="x",main="previous value 7",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)

plotCI(ciriop$V33,xx,li=ciriop$V35,ui=ciriop$V36,xaxt="n",gap=0,ylab="physical integrity",
xlim=c(-1,1),axes=F,err="x",main="previous value 8",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)

dev.off()


### produces Figure 3 

ciriop<-read.csv("ciriopf1sn.csv",header=F,sep="\t")
names(ciriop)
xx<-seq(0:8)

pdf("e:/tex/a/cwhrf1sn.pdf")
postscript("e:/tex/a/cwhrf1sn.eps")
par(mfrow=c(2,3))


plotCI(ciriop$V1,xx,li=ciriop$V3,ui=ciriop$V4,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value < 5",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V5,xx,li=ciriop$V7,ui=ciriop$V8,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 5",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V9,xx,li=ciriop$V11,ui=ciriop$V12,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 6",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V13,xx,li=ciriop$V15,ui=ciriop$V16,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 7",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)


plotCI(ciriop$V17,xx,li=ciriop$V19,ui=ciriop$V20,xaxt="n",gap=0,ylab="physical integrity",
axes=F,xlim=c(-1,1),err="x",main="previous value 8",xlab="change in probability due to
treaty ratification")
 lines(xx*0,xx)
 axis(2, 1:9, 0:8)
  axis(side=1,cex=0.7)

dev.off()

### produces Figure 5
# reads in coefficients and variance-covariance matrix from stata
#matt<-read.dta("e:/info/data/cwhrmat.dta")
matt<-read.dta("cwhrmat.dta")
names(matt)[1:10]
#matt<-matt[,-c(1:3)]
names(matt)
coef<-ncol(matt)/2
matt[1,]
mu<-array(NA,dim=c(1,coef))
mu<-as.numeric(matt[1,1:coef])
sigma<-array(NA,dim=c(coef,coef))
for (i in 1:coef){
sigma[i,]<-as.numeric(matt[i,(coef+1):(2*coef)])
}

str(matt)
coef.sim<-mvrnorm(n=1000,mu,sigma)


rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p<-array(NA,dim=c(9,9,1000))

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 +#cwhr$CompSyst_eu_l1*
 coef.sim[,9]*1 +
 coef.sim[,17]*4740 


p[1,2,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,2,],na.rm=T)

p[2,2,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,2,],na.rm=T)

p[3,2,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,2,],na.rm=T)

p[4,2,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,2,],na.rm=T)

p[5,2,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,2,],na.rm=T)

p[6,2,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,2,],na.rm=T)

p[7,2,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,2,],na.rm=T)

p[8,2,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,2,],na.rm=T)

p[9,2,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,2,],na.rm=T)

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 +#cwhr$CompSyst_eu_l1*
 coef.sim[,10]*1 +
 coef.sim[,18]*4740 


p[1,3,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,3,],na.rm=T)

p[2,3,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,3,],na.rm=T)

p[3,3,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,3,],na.rm=T)

p[4,3,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,3,],na.rm=T)

p[5,3,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,3,],na.rm=T)

p[6,3,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,3,],na.rm=T)

p[7,3,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,3,],na.rm=T)

p[8,3,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,3,],na.rm=T)

p[9,3,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,3,],na.rm=T)

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 +#cwhr$CompSyst_eu_l1*
 coef.sim[,11]*1 +
 coef.sim[,19]*4740 


p[1,4,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,4,],na.rm=T)

p[2,4,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,4,],na.rm=T)

p[3,4,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,4,],na.rm=T)

p[4,4,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,4,],na.rm=T)

p[5,4,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,4,],na.rm=T)

p[6,4,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,4,],na.rm=T)

p[7,4,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,4,],na.rm=T)

p[8,4,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,4,],na.rm=T)

p[9,4,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,4,],na.rm=T)

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 +#cwhr$CompSyst_eu_l1*
 coef.sim[,12]*1 +
 coef.sim[,20]*4740 


p[1,5,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,5,],na.rm=T)

p[2,5,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,5,],na.rm=T)

p[3,5,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,5,],na.rm=T)

p[4,5,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,5,],na.rm=T)

p[5,5,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,5,],na.rm=T)

p[6,5,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,5,],na.rm=T)

p[7,5,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,5,],na.rm=T)

p[8,5,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,5,],na.rm=T)

p[9,5,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,5,],na.rm=T)

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 +#cwhr$CompSyst_eu_l1*
 coef.sim[,13]*1 +
 coef.sim[,21]*4740 


p[1,6,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,6,],na.rm=T)

p[2,6,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,6,],na.rm=T)

p[3,6,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,6,],na.rm=T)

p[4,6,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,6,],na.rm=T)

p[5,6,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,6,],na.rm=T)

p[6,6,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,6,],na.rm=T)

p[7,6,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,6,],na.rm=T)

p[8,6,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,6,],na.rm=T)

p[9,6,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,6,],na.rm=T)

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 +#cwhr$CompSyst_eu_l1*
 coef.sim[,14]*1 +
 coef.sim[,22]*4740 


p[1,7,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,7,],na.rm=T)

p[2,7,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,7,],na.rm=T)

p[3,7,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,7,],na.rm=T)

p[4,7,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,7,],na.rm=T)

p[5,7,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,7,],na.rm=T)

p[6,7,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,7,],na.rm=T)

p[7,7,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,7,],na.rm=T)

p[8,7,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,7,],na.rm=T)

p[9,7,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,7,],na.rm=T)

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 +#cwhr$CompSyst_eu_l1*
 coef.sim[,15]*1 +
 coef.sim[,23]*4740 


p[1,8,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,8,],na.rm=T)

p[2,8,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,8,],na.rm=T)

p[3,8,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,8,],na.rm=T)

p[4,8,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,8,],na.rm=T)

p[5,8,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,8,],na.rm=T)

p[6,8,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,8,],na.rm=T)

p[7,8,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,8,],na.rm=T)

p[8,8,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,8,],na.rm=T)

p[9,8,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,8,],na.rm=T)


rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 0*coef.sim[,4]*0 +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *0 + # *cwhr$RatYears
 coef.sim[,6]*0+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 +#cwhr$CompSyst_eu_l1*
 coef.sim[,16]*1 +
 coef.sim[,24]*4740 


p[1,9,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,9,],na.rm=T)

p[2,9,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,9,],na.rm=T)

p[3,9,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,9,],na.rm=T)

p[4,9,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,9,],na.rm=T)

p[5,9,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,9,],na.rm=T)

p[6,9,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,9,],na.rm=T)

p[7,9,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,9,],na.rm=T)

p[8,9,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,9,],na.rm=T)

p[9,9,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,9,],na.rm=T)

bas<-p


#p1 rat
#p2 rat2
#p3 rat3



rr<-coef.sim[1,1]*cwhr$conf_l1+
 coef.sim[1,2]*cwhr$chga_regime_l1+
 coef.sim[1,2]*cwhr$chga_regime_l1+
 coef.sim[1,3]*cwhr$pennwt_rgdpch_l1 + 
 0*coef.sim[1,4]*(cwhr$Ratification_l1-1) + 
 coef.sim[1,5]    *cwhr$RatYears + 
 coef.sim[1,6]*cwhr$RatYears2 +
 coef.sim[1,7]*cwhr$CompSyst_un_l1+
 coef.sim[1,8]*cwhr$CompSyst_eu_l1+
 coef.sim[1,9]*cwhr$ciri_physint_l11+
 coef.sim[1,10]*cwhr$ciri_physint_l12+
 coef.sim[1,11]*cwhr$ciri_physint_l13+
 coef.sim[1,12]*cwhr$ciri_physint_l14 + 
 0*coef.sim[1,13]* cwhr$ciri_physint_l15+ 
 coef.sim[1,14]    *cwhr$ciri_physint_l16 + 
 coef.sim[1,15]*cwhr$ciri_physint_l17 +
 coef.sim[1,16]**cwhr$ciri_physint_l18

p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *1 + # *cwhr$RatYears
 coef.sim[,6]*1+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *1 + # *cwhr$RatYears
 coef.sim[,6]*1+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *1 + # *cwhr$RatYears
 coef.sim[,60+k]*1 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r1<-p

p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *2 + # *cwhr$RatYears
 coef.sim[,6]*4+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *2 + # *cwhr$RatYears
 coef.sim[,6]*4+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *2 + # *cwhr$RatYears
 coef.sim[,60+k]*4 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r2<-p

p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *3 + # *cwhr$RatYears
 coef.sim[,6]*9+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *3 + # *cwhr$RatYears
 coef.sim[,6]*9+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *3 + # *cwhr$RatYears
 coef.sim[,60+k]*9 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r3<-p

p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *4 + # *cwhr$RatYears
 coef.sim[,6]*16+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *4 + # *cwhr$RatYears
 coef.sim[,6]*16+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *4 + # *cwhr$RatYears
 coef.sim[,60+k]*16 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r4<-p

p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *5 + # *cwhr$RatYears
 coef.sim[,6]*25+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *5 + # *cwhr$RatYears
 coef.sim[,6]*25+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *5 + # *cwhr$RatYears
 coef.sim[,60+k]*25 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r5<-p

p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *6 + # *cwhr$RatYears
 coef.sim[,6]*36+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *6 + # *cwhr$RatYears
 coef.sim[,6]*36+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *6 + # *cwhr$RatYears
 coef.sim[,60+k]*36 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r6<-p

p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *7 + # *cwhr$RatYears
 coef.sim[,6]*49+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *7 + # *cwhr$RatYears
 coef.sim[,6]*49+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *7 + # *cwhr$RatYears
 coef.sim[,60+k]*49 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r7<-p

p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *8 + # *cwhr$RatYears
 coef.sim[,6]*64+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *8 + # *cwhr$RatYears
 coef.sim[,6]*64+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *8 + # *cwhr$RatYears
 coef.sim[,60+k]*64 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r8<-p


p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *9 + # *cwhr$RatYears
 coef.sim[,6]*81+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *9 + # *cwhr$RatYears
 coef.sim[,6]*81+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *9 + # *cwhr$RatYears
 coef.sim[,60+k]*81 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r9<-p


p<-array(NA,dim=c(9,9,1000))

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* +  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *10 + # *cwhr$RatYears
 coef.sim[,6]*100+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*

p[1,1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)

for (k in 1:8){
rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*1+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *10 + # *cwhr$RatYears
 coef.sim[,6]*100+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,44+k]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,52+k]   *10 + # *cwhr$RatYears
 coef.sim[,60+k]*100 #cwhr$RatYears2*0 

p[1,k+1,]<-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[1,1,],na.rm=T)

p[2,k+1,]<-pnorm(coef.sim[,coef-6],mean=rr)-pnorm(coef.sim[,coef-7],mean=rr)
mean(p[2,1,],na.rm=T)

p[3,k+1,]<-pnorm(coef.sim[,coef-5],mean=rr)-pnorm(coef.sim[,coef-6],mean=rr)
mean(p[3,1,],na.rm=T)

p[4,k+1,]<-pnorm(coef.sim[,coef-4],mean=rr)-pnorm(coef.sim[,coef-5],mean=rr)
mean(p[4,1,],na.rm=T)

p[5,k+1,]<-pnorm(coef.sim[,coef-3],mean=rr)-pnorm(coef.sim[,coef-4],mean=rr)
mean(p[5,1,],na.rm=T)

p[6,k+1,]<-pnorm(coef.sim[,coef-2],mean=rr)-pnorm(coef.sim[,coef-3],mean=rr)
mean(p[6,1,],na.rm=T)

p[7,k+1,]<-pnorm(coef.sim[,coef-1],mean=rr)-pnorm(coef.sim[,coef-2],mean=rr)
mean(p[7,1,],na.rm=T)

p[8,k+1,]<-pnorm(coef.sim[,coef],mean=rr)-pnorm(coef.sim[,coef-1],mean=rr)
mean(p[8,1,],na.rm=T)

p[9,k+1,]<-1-pnorm(coef.sim[,coef],mean=rr)
mean(p[9,1,],na.rm=T)
}

r10<-p

rr<-coef.sim[,1]*0+ #cwhr$conf_l1
 coef.sim[,2]*0+ #cwhr$chga_regime_l1*
 coef.sim[,3]*(4730) +  #+0*cwhr$pennwt_rgdpch_l1
 coef.sim[,4]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,5]   *2 + # *cwhr$RatYears
 coef.sim[,6]*4+ #cwhr$RatYears2*0 
 coef.sim[,7]*0+ #*cwhr$CompSyst_un_l1
 coef.sim[,8]*0 #cwhr$CompSyst_eu_l1*
 coef.sim[,45]* 1+  #(cwhr$Ratification_l1-1)*
 coef.sim[,53]   *2 + # *cwhr$RatYears
 coef.sim[,61]*4 #cwhr$RatYears2*0 

p<-array(NA,dim=c(9,9,1000))


# calculates averages and credible intervals

res<-array(NA,dim=c(9,9))
res5<-array(NA,dim=c(9,9))
res95<-array(NA,dim=c(9,9))

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]-bas[,,i]
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}


res01<-res
res015<-res5
res0195<-res95


qq<-p
for (i in 1:1000){
qq[,,i]<-(r1[,,i]%*%r2[,,i])-(bas[,,i]%*%bas[,,i])
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res02<-res
res025<-res5
res0295<-res95

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]%*%r2[,,i]%*%r3[,,i]-bas[,,i]%*%bas[,,i]%*%bas[,,i]
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res03<-res
res035<-res5
res0395<-res95

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]%*%r2[,,i]%*%r3[,,i]%*%r4[,,i]-bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res04<-res
res045<-res5
res0495<-res95

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]%*%r2[,,i]%*%r3[,,i]%*%r4[,,i]%*%r5[,,i]-bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res05<-res
res055<-res5
res0595<-res95

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]%*%r2[,,i]%*%r3[,,i]%*%r4[,,i]%*%r5[,,i]%*%r6[,,i]-bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res06<-res
res065<-res5
res0695<-res95

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]%*%r2[,,i]%*%r3[,,i]%*%r4[,,i]%*%r5[,,i]%*%r6[,,i]%*%r7[,,i]-bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res07<-res
res075<-res5
res0795<-res95

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]%*%r2[,,i]%*%r3[,,i]%*%r4[,,i]%*%r5[,,i]%*%r6[,,i]%*%r7[,,i]%*%r8[,,i]-bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res08<-res
res085<-res5
res0895<-res95

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]%*%r2[,,i]%*%r3[,,i]%*%r4[,,i]%*%r5[,,i]%*%r6[,,i]%*%r7[,,i]%*%r8[,,i]%*%r9[,,i]-bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]
}

for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res09<-res
res095<-res5
res0995<-res95

qq<-p
for (i in 1:1000){
qq[,,i]<-r1[,,i]%*%r2[,,i]%*%r3[,,i]%*%r4[,,i]%*%r5[,,i]%*%r6[,,i]%*%r7[,,i]%*%r8[,,i]%*%r9[,,i]%*%r9[,,i]-bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]%*%bas[,,i]
}


for (i in 1:9){
for (j in 1:9){
res[i,j]<-mean(qq[i,j,],na.rm=T)
res5[i,j]<-quantile(qq[i,j,],probs=seq(0.05,0.05),na.rm=T)
res95[i,j]<-quantile(qq[i,j,],probs=seq(0.95,0.95),na.rm=T)
}}

res10<-res
res105<-res5
res1095<-res95

#pdf("e:/tex/a/cwhrf4n.pdf")
postscript("e:/tex/a/cwhrf4n.eps")
par(mfrow=c(3,4))
par(ps=9)
#par(mfrow=c(1,1))


resp<-rbind(t(res01[,1]),t(res02[,1]),t(res03[,1]),t(res04[,1]),t(res05[,1]),t(res06[,1]),t(res07[,1]),
t(res08[,1]),t(res09[,1]),t(res10[,1]))
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 0",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)


resp<-rbind(t(res01[,2]),t(res02[,2]),t(res03[,2]),t(res04[,2]),t(res05[,2])
,t(res06[,2]),t(res07[,2]),t(res08[,2]),t(res09[,2]),t(res10[,2]))
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 1",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)


resp<-rbind(t(res01[,3]),t(res02[,3]),t(res03[,3]),t(res04[,3]),t(res05[,3])
,t(res06[,3]),t(res07[,3]),t(res08[,3]),t(res09[,3]),t(res10[,3]))
)
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 2",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)

resp<-rbind(t(res01[,4]),t(res02[,4]),t(res03[,4]),t(res04[,4]),t(res05[,4])
,t(res06[,4]),t(res07[,4]),t(res08[,4]),t(res09[,4]),t(res10[,4]))
)
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 3",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)

resp<-rbind(t(res01[,5]),t(res02[,5]),t(res03[,5]),t(res04[,5]),t(res05[,5])
,t(res06[,5]),t(res07[,5]),t(res08[,5]),t(res09[,5]),t(res10[,5]))
)
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 4",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)

resp<-rbind(t(res01[,6]),t(res02[,6]),t(res03[,6]),t(res04[,6]),t(res05[,6])
,t(res06[,6]),t(res07[,6]),t(res08[,6]),t(res09[,6]),t(res10[,6]))
)
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 5",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)

resp<-rbind(t(res01[,7]),t(res02[,7]),t(res03[,7]),t(res04[,7]),t(res05[,7])
,t(res06[,7]),t(res07[,7]),t(res08[,7]),t(res09[,7]),t(res10[,7]))
)
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 6",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)

resp<-rbind(t(res01[,8]),t(res02[,8]),t(res03[,8]),t(res04[,8]),t(res05[,8])
,t(res06[,8]),t(res07[,8]),t(res08[,8]),t(res09[,8]),t(res10[,8]))
)
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 7",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)

resp<-rbind(t(res01[,9]),t(res02[,9]),t(res03[,9]),t(res04[,9]),t(res05[,9])
,t(res06[,9]),t(res07[,9]),t(res08[,9]),t(res09[,9]),t(res10[,9]))
)
xx<-seq(1:nrow(resp))
plot(xx,resp[,1],type="l",ylim=c(0,0),lty=1, main="phys. int. 8",
 xlab="years since ratification",ylab="change in p due to ratific.")
lines(xx,resp[,2],add=T,lty=2)
lines(xx,resp[,3],add=T,lty=3)
lines(xx,resp[,4],add=T,lty=4)
lines(xx,resp[,5],add=T,lty=5)
lines(xx,resp[,6],add=T,lty=6)
lines(xx,resp[,7],add=T,lty=7)
lines(xx,resp[,8],add=T,lty=8)
lines(xx,resp[,9],add=T,lty=9)

plot(0,0,ylim=c(-2.,0.5),xaxt="n",yaxt="n",main="legend",
 xlab="", ylab="",lty=1,type="l")

legend("bottomleft", c("phy.int.0", "phy.int.1", "phy.int.2", "phy.int.3"
, "phy.int.4", "phy.int.5", "phy.int.6", "phy.int.7", "phy.int.8"), 
 lty = c(1, 2, 3,4,5,6,7,8,9)   ,    merge = TRUE, y.intersp=0.75,bty="n")


dev.off()

### produces Figure 5 
 
xxx<-seq(1:11)
#pdf("e:/tex/a/cwhrf3n.pdf")
postscript("e:/tex/a/cwhrf3n.eps")
par(mfrow=c(1,1))

plot(xxx,-0.788+0.215*xxx-0.010*xxx*xxx,type="l",ylim=c(-2.,0.5),
 xlab="years since ratification", ylab="marginal effect on latent variable",lty=1)
lines(xxx,-0.788+0.892+(0.215-0.222)*xxx+(0.011-0.01)*xxx*xxx,type="l",add=T,lty=2)
lines(xxx,-0.788+0.305+(0.215-0.206)*xxx+(0.009-0.01)*xxx*xxx,type="l",add=T,lty=3)
lines(xxx,-0.788+0.194+(0.215-0.141)*xxx+(0.006-0.01)*xxx*xxx,type="l",add=T,lty=4)
lines(xxx,-0.788+0.809+(0.215-0.241)*xxx+(0.010-0.01)*xxx*xxx,type="l",add=T,lty=5)
lines(xxx,-0.788+1.243+(0.215-0.323)*xxx+(0.014-0.01)*xxx*xxx,type="l",add=T,lty=6)
lines(xxx,-0.788+1.350+(0.215-0.278)*xxx+(0.011-0.01)*xxx*xxx,type="l",add=T,lty=7)
lines(xxx,-0.788+0.522+(0.215-0.162)*xxx+(0.007-0.01)*xxx*xxx,type="l",add=T,lty=8)
lines(xxx,-0.788+0.871+(0.215-0.272)*xxx+(0.011-0.01)*xxx*xxx,type="l",add=T,lty=9)
legend("bottomleft", c("phy.int.0", "phy.int.1", "phy.int.2", "phy.int.3"
, "phy.int.4", "phy.int.5", "phy.int.6", "phy.int.7", "phy.int.8"), 
 lty = c(1, 2, 3,4,5,6,7,8,9)   ,    merge = TRUE)

dev.off()

